home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / libs / gle / util / contour / contour.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-29  |  22.6 KB  |  952 lines
  1. /* contour.f -- translated by f2c (version of 26 February 1990  17:38:00).
  2.    You must link the resulting object file with the libraries:
  3.     -lF77 -lI77 -lm -lc   (in that order)
  4. */
  5.  
  6. #include "f2c.h"
  7.  
  8. /* Common Block Declarations */
  9.  
  10. struct {
  11.     real xcur, ycur, xl, yl, cval[10];
  12.     integer clab[10], nch;
  13. } gctcom_;
  14.  
  15. #define gctcom_1 gctcom_
  16.  
  17. /* Table of constant values */
  18.  
  19. static integer c__2 = 2;
  20. static integer c__6 = 6;
  21. static real c_b43 = 0.f;
  22. static real c_b44 = -11.f;
  23. static integer c_n3 = -3;
  24. static real c_b47 = 1.f;
  25. static integer c__1 = 1;
  26. static integer c__51 = 51;
  27. static real c_b55 = 1e6f;
  28. static real c_b56 = .07f;
  29. static integer c__10 = 10;
  30. static integer c__5 = 5;
  31. static real c_b67 = 10.f;
  32. static integer c__999 = 999;
  33. static integer c_n1 = -1;
  34. static integer c__3 = 3;
  35. static real c_b104 = .125f;
  36.  
  37. /* Subroutine */ int fill0_(integer *bitmap, integer *n)
  38. {
  39.     /* Initialized data */
  40.  
  41.     static integer nbpw = 31;
  42.  
  43.     /* System generated locals */
  44.     integer i_1;
  45.  
  46.     /* Builtin functions */
  47.     integer pow_ii(integer *, integer *);
  48.  
  49.     /* Local variables */
  50.     static integer nblw, loop, i;
  51.  
  52.  
  53. /*     fill the first n bits of bitmap with zeroes. */
  54.  
  55.  
  56.     /* Parameter adjustments */
  57.     --bitmap;
  58.  
  59.     /* Function Body */
  60. /*     nbpw is the minimum number of significant bits per word used */
  61. /*     by integer arithmetic.  this is usually one less than the */
  62. /*     actual number of bits per word, but an important exception is */
  63. /*     the cdc-6000 series of machines, where nbpw should be 48. */
  64.  
  65.     loop = *n / nbpw;
  66.     nblw = *n % nbpw;
  67.     if (loop == 0) {
  68.     goto L20;
  69.     }
  70.     i_1 = loop;
  71.     for (i = 1; i <= i_1; ++i) {
  72.     bitmap[i] = 0;
  73. /* L10: */
  74.     }
  75. L20:
  76.     if (nblw != 0) {
  77.     i_1 = nbpw - nblw;
  78.     bitmap[loop + 1] %= pow_ii(&c__2, &i_1);
  79.     }
  80.     return 0;
  81. } /* fill0_ */
  82.  
  83. /* Subroutine */ int mark1_(integer *bitmap, integer *n)
  84. {
  85.     /* Initialized data */
  86.  
  87.     static integer nbpw = 31;
  88.  
  89.     /* System generated locals */
  90.     integer i_1;
  91.  
  92.     /* Builtin functions */
  93.     integer pow_ii(integer *, integer *);
  94.  
  95.     /* Local variables */
  96.     static integer nbit, i, nword;
  97.  
  98.  
  99. /*     put a one in the nth bit of bitmap. */
  100.  
  101.  
  102.     /* Parameter adjustments */
  103.     --bitmap;
  104.  
  105.     /* Function Body */
  106. /*     nbpw is the minimum number of significant bits per word used */
  107. /*     by integer arithmetic.  this is usually one less than the */
  108. /*     actual number of bits per word, but an important exception is */
  109. /*     the cdc-6000 series of machines, where nbpw should be 48. */
  110.  
  111.     nword = (*n - 1) / nbpw;
  112.     nbit = (*n - 1) % nbpw;
  113.     i_1 = nbpw - nbit - 1;
  114.     i = pow_ii(&c__2, &i_1);
  115.     bitmap[nword + 1] += i * (1 - bitmap[nword + 1] / i % 2);
  116.     return 0;
  117. } /* mark1_ */
  118.  
  119. integer iget_(integer *bitmap, integer *n)
  120. {
  121.     /* Initialized data */
  122.  
  123.     static integer nbpw = 31;
  124.  
  125.     /* System generated locals */
  126.     integer ret_val, i_1;
  127.  
  128.     /* Builtin functions */
  129.     integer pow_ii(integer *, integer *);
  130.  
  131.     /* Local variables */
  132.     static integer nbit, nword;
  133.  
  134.  
  135. /*     iget=0 if the nth bit of bitmap is zero, else iget is one. */
  136.  
  137.  
  138.     /* Parameter adjustments */
  139.     --bitmap;
  140.  
  141.     /* Function Body */
  142. /*     nbpw is the minimum number of significant bits per word used */
  143. /*     by integer arithmetic.  this is usually one less than the */
  144. /*     actual number of bits per word, but an important exception is */
  145. /*     the cdc-6000 series of machines, where nbpw should be 48. */
  146.  
  147.     nword = (*n - 1) / nbpw;
  148.     nbit = (*n - 1) % nbpw;
  149.     i_1 = nbpw - nbit - 1;
  150.     ret_val = bitmap[nword + 1] / pow_ii(&c__2, &i_1) % 2;
  151.     return ret_val;
  152. } /* iget_ */
  153.  
  154. /* Subroutine */ int gcontr_(real *z, integer *nrz, integer *nx, integer *ny, 
  155.     real *cv, integer *ncv, real *zmax, integer *bitmap, S_fp draw)
  156. {
  157.     /* Initialized data */
  158.  
  159.     static integer l1[4] = { 0,0,-1,-1 };
  160.     static integer i1[2] = { 1,0 };
  161.     static integer i2[2] = { 1,-1 };
  162.     static integer i3[6] = { 1,0,0,1,1,0 };
  163.  
  164.     /* Format strings */
  165.     static char fmt_100[] = "";
  166.     static char fmt_280[] = "";
  167.  
  168.     /* System generated locals */
  169.     integer z_dim1, z_offset, i_1, i_2, i_3;
  170.     static integer equiv_3[4], equiv_5[2];
  171.     static real equiv_7[2];
  172.  
  173.     /* Builtin functions */
  174.     integer i_sign(integer *, integer *);
  175.  
  176.     /* Local variables */
  177.     static real cval;
  178.     static integer idir;
  179.     extern integer iget_(integer *, integer *);
  180.     static real dmax_;
  181. #define imin (equiv_3 + 2)
  182. #define imax (equiv_3)
  183. #define jmax (equiv_3 + 1)
  184. #define jmin (equiv_3 + 3)
  185.     static integer icur, jcur, jump;
  186.     static real xint[4];
  187.     extern /* Subroutine */ int fill0_(integer *, integer *), mark1_(integer *
  188.         , integer *);
  189. #define i (equiv_5)
  190. #define j (equiv_5 + 1)
  191.     static integer k, l, iedge, iflag;
  192. #define x (equiv_7)
  193. #define y (equiv_7 + 1)
  194.     static integer ibkey;
  195. #define l2 (equiv_3)
  196.     static real z1, z2;
  197.     static integer ii;
  198. #define ij (equiv_5)
  199.     static integer jj, ni, ks, ix;
  200. #define xy (equiv_7)
  201.     static real zz;
  202.     static integer nxidir, icv;
  203.  
  204.     /* Assigned format variables */
  205.     char *jump_fmt;
  206.  
  207.  
  208. /*     this subroutine draws a contour through equal values of an array. 
  209. */
  210.  
  211. /*     *****     formal arguments     *********************************** 
  212. */
  213.  
  214. /*     z is the array for which contours are to be drawn.  the elements */
  215.  
  216. /*     of z are assumed to lie upon the nodes of a topologically */
  217. /*     rectangular coordinate system - e.g. cartesian, polar (except */
  218. /*     the origin), etc. */
  219.  
  220. /*     nrz is the number of rows declared for z in the calling program. */
  221.  
  222.  
  223. /*     nx is the limit for the first subscript of z. */
  224.  
  225. /*     ny is the limit for the second subscript of z. */
  226.  
  227. /*     cv are the values of the contours to be drawn. */
  228.  
  229. /*     ncv is the number of contour values in cv. */
  230.  
  231. /*     zmax is the maximum value of z for consideration.  a value of */
  232. /*     z(i,j) greater than zmax is a signal that that point and the */
  233. /*     grid line segments radiating from that point to it's neighbors */
  234. /*     are to be excluded from contouring. */
  235.  
  236. /*     bitmap is a work area large enough to hold 2*nx*ny*ncv bits.  it */
  237.  
  238. /*     is accessed by low-level routines, which are described below. */
  239. /*     let j be the number of useful bits in each word of bitmap, */
  240. /*     as determined by the user machine and implementation of */
  241. /*     the bitmap manipulation subprograms described below.  then */
  242. /*     the number of words required for the bitmap is the floor of */
  243. /*         (2*nx*ny*ncv+j-1)/j. */
  244.  
  245. /*     draw is a user-provided subroutine used to draw contours. */
  246. /*     the calling sequence for draw is: */
  247.  
  248. /*         call draw (x,y,iflag) */
  249. /*         let nx = integer part of x, fx = fractional part of x. */
  250. /*         then x should be interpreted such that increases in nx */
  251. /*         correspond to increases in the first subscript of z, and */
  252. /*         fx is the fractional distance from the abscissa corresponding 
  253. */
  254. /*         to nx to the abscissa corresponding to nx+1, */
  255. /*         and y should be interpreted similarly for the second */
  256. /*         subscript of z. */
  257. /*         the low-order digit of iflag will have one of the values: */
  258. /*             1 - continue a contour, */
  259. /*             2 - start a contour at a boundary, */
  260. /*             3 - start a contour not at a boundary, */
  261. /*             4 - finish a contour at a boundary, */
  262. /*             5 - finish a closed contour (not at a boundary). */
  263. /*                 note that requests 1, 4 and 5 are for pen-down */
  264. /*                 moves, and that requests 2 and 3 are for pen-up */
  265. /*                 moves. */
  266. /*             6 - set x and y to the approximate 'pen' position, using */
  267.  
  268. /*                 the notation discussed above.  this call may be */
  269. /*                 ignored, the result being that the 'pen' position */
  270. /*                 is taken to correspond to z(1,1). */
  271. /*         iflag/10 is the contour number. */
  272.  
  273. /*     *****     external subprograms     ******************************* 
  274. */
  275.  
  276. /*     draw is the user-supplied line drawing subprogram described above. 
  277. */
  278. /*     draw may be sensitive to the host computer and to the plot device. 
  279. */
  280. /*     fill0 is used to fill a bitmap with zeroes.  call fill0 (bitmap,n) 
  281. */
  282. /*     fills the first n bits of bitmap with zeroes. */
  283. /*     mark1 is used to place a 1 in a specific bit of the bitmap. */
  284. /*     call mark1 (bitmap,n) puts a 1 in the nth bit of the bitmap. */
  285. /*     iget is used to determine the setting of a particular bit in the */
  286.  
  287. /*     bitmap.  i=iget(bitmap,n) sets i to zero if the nth bit of the */
  288. /*     bitmap is zero, and sets i to one if the nth bit is one. */
  289. /*     fill0, mark1 and iget are machine sensitive. */
  290.  
  291. /*     ****************************************************************** 
  292. */
  293.  
  294.  
  295. /*     l1 and l2 contain limits used during the spiral search for the */
  296. /*     beginning of a contour. */
  297. /*     ij stores subcripts used during the spiral search. */
  298.  
  299.  
  300. /*     i1, i2 and i3 are used for subscript computations during the */
  301. /*     examination of lines from z(i,j) to it's neighbors. */
  302.  
  303.  
  304. /*     xint is used to mark intersections of the contour under */
  305. /*     consideration with the edges of the cell being examined. */
  306.  
  307.  
  308. /*     xy is used to compute coordinates for the draw subroutine. */
  309.  
  310.  
  311.     /* Parameter adjustments */
  312.     --bitmap;
  313.     --cv;
  314.     z_dim1 = *nrz;
  315.     z_offset = z_dim1 + 1;
  316.     z -= z_offset;
  317.  
  318.     /* Function Body */
  319.  
  320.     l1[0] = *nx;
  321.     l1[1] = *ny;
  322.     dmax_ = *zmax;
  323.  
  324. /*     set the current pen position.  the default position corresponds */
  325. /*     to z(1,1). */
  326.  
  327.     *x = 1.f;
  328.     *y = 1.f;
  329.     (*draw)(x, y, &c__6);
  330. /* Computing MAX */
  331. /* Computing MIN */
  332.     i_3 = (integer) (*x);
  333.     i_1 = 1, i_2 = min(i_3,*nx);
  334.     icur = max(i_1,i_2);
  335. /* Computing MAX */
  336. /* Computing MIN */
  337.     i_3 = (integer) (*y);
  338.     i_1 = 1, i_2 = min(i_3,*ny);
  339.     jcur = max(i_1,i_2);
  340.  
  341. /*     clear the bitmap */
  342.  
  343.     i_1 = (*nx << 1) * *ny * *ncv;
  344.     fill0_(&bitmap[1], &i_1);
  345.  
  346. /*     search along a rectangular spiral path for a line segment having */
  347.  
  348. /*     the following properties: */
  349. /*          1.  the end points are not excluded, */
  350. /*          2.  no mark has been recorded for the segment, */
  351. /*          3.  the values of z at the ends of the segment are such that 
  352. */
  353. /*              one z is less than the current contour value, and the */
  354. /*              other is greater than or equal to the current contour */
  355. /*              value. */
  356.  
  357. /*     search all boundaries first, then search interior line segments. */
  358.  
  359. /*     note that the interior line segments near excluded points may be */
  360.  
  361. /*     boundaries. */
  362.  
  363.     ibkey = 0;
  364. L10:
  365.     *i = icur;
  366.     *j = jcur;
  367. L20:
  368.     *imax = *i;
  369.     *imin = -(*i);
  370.     *jmax = *j;
  371.     *jmin = -(*j);
  372.     idir = 0;
  373. /*     direction zero is +i, 1 is +j, 2 is -i, 3 is -j. */
  374. L30:
  375.     nxidir = idir + 1;
  376.     k = nxidir;
  377.     if (nxidir > 3) {
  378.     nxidir = 0;
  379.     }
  380. L40:
  381.     *i = abs(*i);
  382.     *j = abs(*j);
  383.     if (z[*i + *j * z_dim1] > dmax_) {
  384.     goto L140;
  385.     }
  386.     l = 1;
  387. /*     l=1 means horizontal line, l=2 means vertical line. */
  388. L50:
  389.     if (ij[l - 1] >= l1[l - 1]) {
  390.     goto L130;
  391.     }
  392.     ii = *i + i1[l - 1];
  393.     jj = *j + i1[3 - l - 1];
  394.     if (z[ii + jj * z_dim1] > dmax_) {
  395.     goto L130;
  396.     }
  397.     jump = 0;
  398.     jump_fmt = fmt_100;
  399. /*     the next 15 statements (or so) detect boundaries. */
  400. L60:
  401.     ix = 1;
  402.     if (ij[3 - l - 1] == 1) {
  403.     goto L80;
  404.     }
  405.     ii = *i - i1[3 - l - 1];
  406.     jj = *j - i1[l - 1];
  407.     if (z[ii + jj * z_dim1] > dmax_) {
  408.     goto L70;
  409.     }
  410.     ii = *i + i2[l - 1];
  411.     jj = *j + i2[3 - l - 1];
  412.     if (z[ii + jj * z_dim1] < dmax_) {
  413.     ix = 0;
  414.     }
  415. L70:
  416.     if (ij[3 - l - 1] >= l1[3 - l - 1]) {
  417.     goto L90;
  418.     }
  419. L80:
  420.     ii = *i + i1[3 - l - 1];
  421.     jj = *j + i1[l - 1];
  422.     if (z[ii + jj * z_dim1] > dmax_) {
  423.     goto L90;
  424.     }
  425.     if (z[*i + 1 + (*j + 1) * z_dim1] < dmax_) {
  426.     switch (jump) {
  427.         case 0: goto L100;
  428.         case 1: goto L280;
  429.     }
  430.     }
  431. L90:
  432.     ix += 2;
  433.     switch (jump) {
  434.     case 0: goto L100;
  435.     case 1: goto L280;
  436.     }
  437. L100:
  438.     if (ix == 3) {
  439.     goto L130;
  440.     }
  441.     if (ix + ibkey == 0) {
  442.     goto L130;
  443.     }
  444. /*     now determine whether the line segment is crossed by the contour. 
  445. */
  446.     ii = *i + i1[l - 1];
  447.     jj = *j + i1[3 - l - 1];
  448.     z1 = z[*i + *j * z_dim1];
  449.     z2 = z[ii + jj * z_dim1];
  450.     i_1 = *ncv;
  451.     for (icv = 1; icv <= i_1; ++icv) {
  452.     i_2 = (*nx * (*ny * (icv - 1) + *j - 1) + *i - 1 << 1) + l;
  453.     if (iget_(&bitmap[1], &i_2) != 0) {
  454.         goto L120;
  455.     }
  456.     if (cv[icv] <= dmin(z1,z2)) {
  457.         goto L110;
  458.     }
  459.     if (cv[icv] <= dmax(z1,z2)) {
  460.         goto L190;
  461.     }
  462. L110:
  463.     i_2 = (*nx * (*ny * (icv - 1) + *j - 1) + *i - 1 << 1) + l;
  464.     mark1_(&bitmap[1], &i_2);
  465. L120:
  466.     ;}
  467. L130:
  468.     ++l;
  469.     if (l <= 2) {
  470.     goto L50;
  471.     }
  472. L140:
  473.     l = idir % 2 + 1;
  474.     ij[l - 1] = i_sign(&ij[l - 1], &l1[k - 1]);
  475.  
  476. /*     lines from z(i,j) to z(i+1,j) and z(i,j+1) are not satisfactory. */
  477.  
  478. /*     continue the spiral. */
  479.  
  480. L150:
  481.     if (ij[l - 1] >= l1[k - 1]) {
  482.     goto L170;
  483.     }
  484.     ++ij[l - 1];
  485.     if (ij[l - 1] > l2[k - 1]) {
  486.     goto L160;
  487.     }
  488.     goto L40;
  489. L160:
  490.     l2[k - 1] = ij[l - 1];
  491.     idir = nxidir;
  492.     goto L30;
  493. L170:
  494.     if (idir == nxidir) {
  495.     goto L180;
  496.     }
  497.     ++nxidir;
  498.     ij[l - 1] = l1[k - 1];
  499.     k = nxidir;
  500.     l = 3 - l;
  501.     ij[l - 1] = l2[k - 1];
  502.     if (nxidir > 3) {
  503.     nxidir = 0;
  504.     }
  505.     goto L150;
  506. L180:
  507.     if (ibkey != 0) {
  508.     return 0;
  509.     }
  510.     ibkey = 1;
  511.     goto L10;
  512.  
  513. /*     an acceptable line segment has been found. */
  514. /*     follow the contour until it either hits a boundary or closes. */
  515.  
  516. L190:
  517.     iedge = l;
  518.     cval = cv[icv];
  519.     if (ix != 1) {
  520.     iedge += 2;
  521.     }
  522.     iflag = ibkey + 2;
  523.     xint[iedge - 1] = (cval - z1) / (z2 - z1);
  524. L200:
  525.     xy[l - 1] = (real) ij[l - 1] + xint[iedge - 1];
  526.     xy[3 - l - 1] = (real) ij[3 - l - 1];
  527.     i_1 = (*nx * (*ny * (icv - 1) + *j - 1) + *i - 1 << 1) + l;
  528.     mark1_(&bitmap[1], &i_1);
  529.     i_1 = iflag + icv * 10;
  530.     (*draw)(x, y, &i_1);
  531.     if (iflag < 4) {
  532.     goto L210;
  533.     }
  534.     icur = *i;
  535.     jcur = *j;
  536.     goto L20;
  537.  
  538. /*     continue a contour.  the edges are numbered clockwise with */
  539. /*     the bottom edge being edge number one. */
  540.  
  541. L210:
  542.     ni = 1;
  543.     if (iedge < 3) {
  544.     goto L220;
  545.     }
  546.     *i -= i3[iedge - 1];
  547.     *j -= i3[iedge + 1];
  548. L220:
  549.     for (k = 1; k <= 4; ++k) {
  550.     if (k == iedge) {
  551.         goto L250;
  552.     }
  553.     ii = *i + i3[k - 1];
  554.     jj = *j + i3[k];
  555.     z1 = z[ii + jj * z_dim1];
  556.     ii = *i + i3[k];
  557.     jj = *j + i3[k + 1];
  558.     z2 = z[ii + jj * z_dim1];
  559.     if (cval <= dmin(z1,z2)) {
  560.         goto L250;
  561.     }
  562.     if (cval > dmax(z1,z2)) {
  563.         goto L250;
  564.     }
  565.     if (k == 1) {
  566.         goto L230;
  567.     }
  568.     if (k != 4) {
  569.         goto L240;
  570.     }
  571. L230:
  572.     zz = z1;
  573.     z1 = z2;
  574.     z2 = zz;
  575. L240:
  576.     xint[k - 1] = (cval - z1) / (z2 - z1);
  577.     ++ni;
  578.     ks = k;
  579. L250:
  580.     ;}
  581.     if (ni == 2) {
  582.     goto L260;
  583.     }
  584.  
  585. /*     the contour crosses all four edges of the cell being examined. */
  586. /*     choose the lines top-to-left and bottom-to-right if the */
  587. /*     interpolation point on the top edge is less than the interpolation 
  588. */
  589. /*     point on the bottom edge.  otherwise, choose the other pair.  this 
  590. */
  591. /*     method produces the same results if the axes are reversed.  the */
  592. /*     contour may close at any edge, but must not cross itself inside */
  593. /*     any cell. */
  594.  
  595.     ks = 5 - iedge;
  596.     if (xint[2] < xint[0]) {
  597.     goto L260;
  598.     }
  599.     ks = 3 - iedge;
  600.     if (ks <= 0) {
  601.     ks += 4;
  602.     }
  603.  
  604. /*     determine whether the contour will close or run into a boundary */
  605. /*     at edge ks of the current cell. */
  606.  
  607. L260:
  608.     l = ks;
  609.     iflag = 1;
  610.     jump = 1;
  611.     jump_fmt = fmt_280;
  612.     if (ks < 3) {
  613.     goto L270;
  614.     }
  615.     *i += i3[ks - 1];
  616.     *j += i3[ks + 1];
  617.     l = ks - 2;
  618. L270:
  619.     i_1 = (*nx * (*ny * (icv - 1) + *j - 1) + *i - 1 << 1) + l;
  620.     if (iget_(&bitmap[1], &i_1) == 0) {
  621.     goto L60;
  622.     }
  623.     iflag = 5;
  624.     goto L290;
  625. L280:
  626.     if (ix != 0) {
  627.     iflag = 4;
  628.     }
  629. L290:
  630.     iedge = ks + 2;
  631.     if (iedge > 4) {
  632.     iedge += -4;
  633.     }
  634.     xint[iedge - 1] = xint[ks - 1];
  635.     goto L200;
  636.  
  637. } /* gcontr_ */
  638.  
  639. #undef xy
  640. #undef ij
  641. #undef l2
  642. #undef y
  643. #undef x
  644. #undef j
  645. #undef i
  646. #undef jmin
  647. #undef jmax
  648. #undef imax
  649. #undef imin
  650.  
  651.  
  652. /*      dimension z(51,51), c(10), work(1680) */
  653. /*     dimension of work is large enough to contain */
  654. /*     2*(dimension of c)*(total dimension of z) useful bits.  see the */
  655. /*     bitmap routines accessed by gcontr. */
  656. /*      real mu */
  657. /*      external draw */
  658. /*      common /cur/ xcur, ycur */
  659. /*      data c(1), c(2), c(3), c(4), c(5) /3.05,3.2,3.5,3.50135,3.6/ */
  660. /*      data c(6), c(7), c(8), c(9), c(10) /3.766413,4.0,4.130149,5.0, */
  661. /*     *  10.0/ */
  662. /*      data nx /51/, ny /51/, nf /10/ */
  663. /*      data xmin /-2.0/, xmax /2.0/, ymin /-2.0/, ymax /2.0/, mu /0.3/ */
  664. /*      dx = (xmax-xmin)/float(nx-1) */
  665. /*      dy = (ymax-ymin)/float(ny-1) */
  666. /*      xcur = 1.0 */
  667. /*      ycur = 1.0 */
  668. /*      if (mod(nx,2).ne.0) ycur = float(ny) */
  669. /*      if (mod(ny,2).ne.0) xcur = float(nx) */
  670. /*      x = xmin - dx */
  671. /*      do 20 i=1,nx */
  672. /*        y = ymin - dy */
  673. /*        x = x + dx */
  674. /*        do 10 j=1,ny */
  675. /*          y = y + dy */
  676. /*          z(i,j) = (1.0-mu)*(2.0/sqrt((x-mu)**2+y**2)+(x-mu)**2+y**2) */
  677. /*     *      + mu*(2.0/sqrt((x+1.0-mu)**2+y**2)+(x+1.0-mu)**2+y**2) */
  678. /*   10   continue */
  679. /*   20 continue */
  680. /*      call gcontr(z, 51, nx, ny, c, nf, 1.e6, work, draw) */
  681. /*      stop */
  682. /*      end */
  683.  
  684. /*      subroutine draw(x, y, iflag) */
  685.  
  686. /*     do output for gcontr. */
  687.  
  688. /*      integer print */
  689. /*      common /cur/ xcur, ycur */
  690. /*      data print /6/ */
  691. /*     print is the system printer fortran i/o unit number. */
  692. /*      icont = iflag/10 */
  693. /*      jump = mod(iflag,10) */
  694. /*      go to (10, 20, 30, 40, 50, 60), jump */
  695. /*   10 write (print,99999) icont, x, y */
  696. /*      go to 70 */
  697. /*   20 write (print,99998) icont, x, y */
  698. /*      go to 70 */
  699. /*   30 write (print,99997) icont, x, y */
  700. /*      go to 70 */
  701. /*   40 write (print,99996) icont, x, y */
  702. /*      go to 70 */
  703. /*   50 write (print,99995) icont, x, y */
  704. /*      go to 70 */
  705. /*   60 write (print,99994) */
  706. /*      x = xcur */
  707. /*      y = ycur */
  708. /*   70 return */
  709. /* 99999 format (17h continue contour, i3, 3h to, 1p2e14.7) */
  710. /* 99998 format (14h start contour, i3, 19h on the boundary at, 1p2e14.7) */
  711. /* 99997 format (14h start contour, i3, 19h in the interior at, 1p2e14.7) */
  712. /* 99996 format (15h finish contour, i3, 19h on the boundary at, 1p2e14.7) */
  713. /* 99995 format (15h finish contour, i3, 19h in the interior at, 1p2e14.7) */
  714. /* 99994 format (33h request for current pen position) */
  715. /*      end */
  716. /* Subroutine */ int cgrid_(integer *nopt, integer *nx, real *sx, real *xs, 
  717.     real *xf, integer *ny, real *sy, real *ys, real *yf)
  718. {
  719.     /* System generated locals */
  720.     integer i_1;
  721.     real r_1, r_2;
  722.  
  723.     /* Builtin functions */
  724.     double r_sign(real *, real *);
  725.  
  726.     /* Local variables */
  727.     static real xinc, yinc, xmin, ymin, xmax, ymax;
  728.     extern /* Subroutine */ int plot_(real *, real *, integer *);
  729.     static integer i, j, n;
  730.     static real xlgth, ylgth, x1, x2, y2, y1, xx;
  731.  
  732.  
  733. /* subroutine which draws a frame around the plot and draws */
  734. /* either tick marks or grid lines. */
  735.  
  736. /* parameters:  nopt -- =0, draw ticks only */
  737. /*                      =1, draw grid lines */
  738. /*                      =2, draw grid lines to edge of frame. */
  739. /*              nx -- number of intervals in x direction */
  740. /*              sx -- spacing in inches between tick marks or grid lines 
  741. */
  742. /*                    along the x axis */
  743. /*              xs -- location of first tick or grid line on x axis */
  744. /*              xf -- location of right edge of frame */
  745. /*              ny -- number of intervals in y direction */
  746. /*              sy -- spacing in inches between tick marks or grid lines 
  747. */
  748. /*                    along the y axis */
  749. /*              ys -- location of first tick or grid line on y axis */
  750. /*              yf -- location of top edge of frame */
  751. /* assumptions: nx, sx, ny, sy all positive. */
  752. /*              the lower left-hand corner of the frame is drawn at (0,0) 
  753. */
  754. /*              if xs<0, use 0; if ys<0, use 0 */
  755. /*              if xf<=0, use nx*sx; if yf<=0, use ny*sy. */
  756.  
  757.     xinc = *sx;
  758.     yinc = *sy;
  759.     xlgth = (real) (*nx) * *sx;
  760.     ylgth = (real) (*ny) * *sy;
  761.     xmin = dmax(*xs,0.f);
  762.     ymin = dmax(*ys,0.f);
  763. /* Computing MAX */
  764.     r_1 = *xf, r_2 = xlgth + xmin;
  765.     xmax = dmax(r_1,r_2);
  766. /* Computing MAX */
  767.     r_1 = *yf, r_2 = ylgth + ymin;
  768.     ymax = dmax(r_1,r_2);
  769.  
  770. /*     draw frame. */
  771.  
  772.     plot_(&c_b43, &c_b43, &c__3);
  773.     plot_(&xmax, &c_b43, &c__2);
  774.     plot_(&xmax, &ymax, &c__2);
  775.     plot_(&c_b43, &ymax, &c__2);
  776.     plot_(&c_b43, &c_b43, &c__2);
  777.     if (*nopt != 0) {
  778.     goto L130;
  779.     }
  780.  
  781. /*     draw tick marks. */
  782.  
  783.     for (j = 1; j <= 4; ++j) {
  784.     switch (j) {
  785.         case 1:  goto L10;
  786.         case 2:  goto L50;
  787.         case 3:  goto L20;
  788.         case 4:  goto L40;
  789.     }
  790. L10:
  791.     x2 = 0.f;
  792.     if (xmin != 0.f) {
  793.         x2 = xmin - *sx;
  794.     }
  795.     y2 = 0.f;
  796.     goto L30;
  797. L20:
  798.     xinc = -(doublereal)(*sx);
  799.     x2 = xmin + xlgth + *sx;
  800.     if (xmax == xmin + xlgth) {
  801.         x2 = xmax;
  802.     }
  803.     y2 = ymax;
  804. L30:
  805.     y1 = y2;
  806.     y2 += r_sign(&c_b104, &xinc);
  807.     n = *nx;
  808.     if ((r_1 = xmax - xmin - xlgth, dabs(r_1)) + dabs(xmin) != 0.f) {
  809.         goto L70;
  810.     } else {
  811.         goto L80;
  812.     }
  813. L40:
  814.     yinc = -(doublereal)(*sy);
  815.     y2 = ymin + ylgth + *sy;
  816.     if (ymax == ymin + ylgth) {
  817.         y2 = ymax;
  818.     }
  819.     x2 = 0.f;
  820.     goto L60;
  821. L50:
  822.     y2 = 0.f;
  823.     if (ymin != 0.f) {
  824.         y2 = ymin - *sy;
  825.     }
  826.     x2 = xmax;
  827. L60:
  828.     x1 = x2;
  829.     n = *ny;
  830.     x2 -= r_sign(&c_b104, &yinc);
  831.     if ((r_1 = ymax - ymin - ylgth, dabs(r_1)) + dabs(ymin) != 0.f) {
  832.         goto L70;
  833.     } else {
  834.         goto L80;
  835.     }
  836. L70:
  837.     ++n;
  838. L80:
  839.     i_1 = n;
  840.     for (i = 1; i <= i_1; ++i) {
  841.         if (j % 2 == 0) {
  842.         goto L90;
  843.         }
  844.         x2 += xinc;
  845.         x1 = x2;
  846.         goto L100;
  847. L90:
  848.         y2 += yinc;
  849.         y1 = y2;
  850. L100:
  851.         plot_(&x1, &y1, &c__3);
  852.         plot_(&x2, &y2, &c__2);
  853. /* L110: */
  854.     }
  855. /* L120: */
  856.     }
  857.     goto L240;
  858.  
  859. /*     draw grid lines */
  860.  
  861. L130:
  862.     x1 = xmin;
  863.     x2 = xmin + xlgth;
  864.     if (*nopt != 2) {
  865.     goto L140;
  866.     }
  867.     x1 = 0.f;
  868.     x2 = xmax;
  869. L140:
  870.     y1 = ymin - *sy;
  871.     n = *ny + 1;
  872.     if (ymax == ymin + ylgth) {
  873.     --n;
  874.     }
  875.     if (ymin != 0.f) {
  876.     goto L150;
  877.     }
  878.     y1 = 0.f;
  879.     --n;
  880. L150:
  881.     if (n <= 0) {
  882.     goto L170;
  883.     }
  884.     j = 1;
  885.     i_1 = n;
  886.     for (i = 1; i <= i_1; ++i) {
  887.     j = -j;
  888.     y1 += *sy;
  889.     plot_(&x1, &y1, &c__3);
  890.     plot_(&x2, &y1, &c__2);
  891.     xx = x1;
  892.     x1 = x2;
  893.     x2 = xx;
  894. /* L160: */
  895.     }
  896. L170:
  897.     y1 = ymin + ylgth;
  898.     y2 = ymin;
  899.     if (*nopt != 2) {
  900.     goto L180;
  901.     }
  902.     y1 = ymax;
  903.     y2 = 0.f;
  904. L180:
  905.     n = *nx + 1;
  906.     if (j < 0) {
  907.     goto L200;
  908.     }
  909.     x1 = xmin - *sx;
  910.     if (xmax == xmin + xlgth) {
  911.     --n;
  912.     }
  913.     if (xmin != 0.f) {
  914.     goto L190;
  915.     }
  916.     x1 = 0.f;
  917.     --n;
  918. L190:
  919.     if (n <= 0) {
  920.     goto L240;
  921.     }
  922.     xinc = *sx;
  923.     goto L220;
  924. L200:
  925.     x1 = xmin + xlgth + *sx;
  926.     if (xmin == 0.f) {
  927.     --n;
  928.     }
  929.     if (xmax != xlgth + xmin) {
  930.     goto L210;
  931.     }
  932.     --n;
  933.     x1 = xmax;
  934. L210:
  935.     xinc = -(doublereal)(*sx);
  936. L220:
  937.     i_1 = n;
  938.     for (i = 1; i <= i_1; ++i) {
  939.     x1 += xinc;
  940.     plot_(&x1, &y1, &c__3);
  941.     plot_(&x1, &y2, &c__2);
  942.     xx = y1;
  943.     y1 = y2;
  944.     y2 = xx;
  945. /* L230: */
  946.     }
  947. L240:
  948.     return 0;
  949.  
  950. } /* cgrid_ */
  951.  
  952.